Rugged Fitness Landscapes of Kauffman Models with a Scale-Free Network 
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We study the nature of fitness landscapes of 'quenched' Kauffman's Boolean model with a scale- 
free network. We have numerically calculated the rugged fitness landscapes, the distributions, its 
tails, and the correlation between the fitness of local optima and their Hamming distance from 
the highest optimum found, respectively. We have found that (a) there is an interesting difference 
between the random and the scale-free networks such that the statistics of the rugged fitness land- 
scapes is Gaussian for the random network while it is non-Gaussian with a tail for the scale-free 
network; (b) as the average degree (k) increases, there is a phase transition at the critical value of 
(k) — (k) c — 2, below which there is a global order and above which the order goes away. 
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I. INTRODUCTION 

The origin of life is one of the most important unsolved 
problems in science lj. To answer the quest, the self- 
organization of matterQ and the emergence of order[^ 
have been regarded as the key ideas. To investigate such 
ideas, as early as in 1969 Kauffman has introduced the so- 
called NK Kauffman model - a random Boolean network 
model based upon the random network theory 4J . This 
model has been a prototype model and studied by many 
authors for a long time in understanding complex systems 
such as metabolic stability and epigenesis, genetic regula- 
tory networks, and transcriptional networks 3| as well as 
general Boolean networks^ 0, 0, 0, neural networks^ 
and spin glasses jToj . 

At nearly the end of 1990's new kinds of networks, 
called scale-free networks, have been discovered from 

studying the growth of the internet geometry and 

topology [Til Il2t fl3| . After the discovery, scientists have 
known that many systems such as those which were orig- 
inally studied by Kauffman as well as other various sys- 
tems such as internet topology, human sexual relation- 
ship, scientific collaboration, economical network, etc. 
belong to the category of the scale-free networks. There- 
fore, it is very interesting for us to know what will happen 
when we apply the concepts of scale-free networks to the 
Kauffman's Boolean network models. 

Recentl y, t here have just started some studies in this 
direction |l3, EE EE E3, EE O E3, EH. There, the 
Boolean dynamics of the Kauffman model with a scale- 
free network has been intensively studied. We would like 
to shortly call this the scale-free (SF) Kauffman model. 
Hence, there are many interesting problems that are nec- 
essary to be considered. 

Thus, in this paper we would like to study the struc- 
ture of rugged fitness landscapes of the SF Kauffman 
model. We would like to know what is the difference in 
fitness landscapes between the NK and the SF Kauff- 
man models. 



The organization of the paper is the following: In 
SecJI, we summarize the formalism on Boolean dynam- 
ics of the NK Kauffman model. In Sec. Ill, we intro- 
duce the formalism to calculate the statistics and the 
rugged fitness landscapes of NK Kauffman model. And 
our scheme to obtain a scale-free network is given. In 
Sec. IV, we show the numerical results of the rugged fit- 
ness landscapes, of the histograms and its tails, and of the 
correlation between the fitness of local optima and their 
Hamming distance from the highest optimum found, for 
both 'quenched' NK and SF models, respectively. In 
Sec.V, our conclusion will be given. 



II. NK KAUFFMAN MODEL 

In the NK model (the SF model is described below), 
we assume that the total number of nodes (vertices) N 
and the degree (i.e., the number of inputs) of the i-th 
node ki in the network are fixed such that all fcj = K. 
Therefore, the resultant graph is a directed random net- 
work, where each link has its own direction as represented 
by an arrow on the link. This gives us in general an 
asymmetric adjacency matrix of network theory. Since 
there are K inputs to each node, 2 K Boolean spin con- 
figurations can be defined on each node; the number 2 K 
certainly becomes very large as K becomes a large num- 
ber. 

We then assume that Boolean functions are randomly 
chosen on each node from the 2 2 possibilities. Locally 
this can be represented by 
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for i = 1, . . . , N, where ct, £ Z% = {0, 1} is the binary 
state and Bi £ Zi is a Boolean function at the «th node, 
randomly chosen from 2 2 * Boolean functions with the 
probability p (or 1 — p) to take 1 (or 0). 

If we fixed the set of the randomly chosen Boolean 
functions {Bi,i = 1, ••■ ,N} in the course of the time 



2 



development, then this model is called the quenched 
model^. On the other hand, if we change the set each 
time, then this model is called the annealed model 
If we study the dynamics of the states in the system tak- 
ing care of Eq.(l), then we are able to obtain the cyclic 
structures of the states such as the length of the cycle, the 
transient time and the basin sizes, etc. These are usually 
calculated numerically, since it is extremely difficult to 
do the calculations analytically [!| 

However, in the annealed models |j.|tjj. it has been in- 
vestigated analytically that as the degree of nodes K is 
increasing, there exists a kind of phase transition of net- 
work at the critical degree K c = \/[2p(l — p)\ and if we 
conversely solve it for p then we obtain the critical prob- 
ability p c = (l± y/i - 2/K)/2. 

III. FITNESS LANDSCAPE MODEL 

Let us now study statistics in the structure of the fit- 
ness landscapes of the NK and SF Kauffman models. 
The fitness landscapes are calculated as follows 0: (i) 
Generate a network with N nodes and the degree fej of 
the i-th node. The ki links are inputs that are directed 
to the i-th node, (ii) Define the local fitness at the i-th 
node by 

Wi = fi{<Ji)(Ti^- ■ ■ ,a lk J (2) 

for % — 1,...,N, where Wi takes one of 2 i+1 real 
numbers u>i which are randomly taken from the inter- 
val [0,1]. This provides us a table for each node (TA- 
BLET), (iii) Define an TV-component initial state A, say 
= (0, 1, 1, . . . , 0). (iv) Investigate the input states 
on ki links for the i-th node. And adjusting the states 
in the entries with the table, choose the fitness Wi from 
the 2 ki+1 values of Wi, u>2, ■ ■ • , u> 2 fci+i. Then, define the 
fitness Wa for the state A by 

1 N 

i=l 

(v) Each state forms a vertex of the iV-dimensional hy- 
percube so that there are totally 2 N vertices, each of 
which has its N neighbors such as ^b = (1, 1,1, .. .,0) 
(i.e., one-mutant variants, denoted by n omv ). Then, cal- 
culate the fitnesses Wb for these neighbor states in the 
same way. (vi) Compare the fitness value Wa of the 
A state with those of the neighbor states such as Wb, 
successively. Here the Hamming distances between the 
state and its neighbors are all 1. If Wa > Wb for all 
neighbors B's, then the fitness Wa for the A state is a 
local optimum. And if we meet a neighbor B such that 
Wa<W b , then write B — A' and W B = W A >- Redo the 
same procedure with all neighbors to obtain Wa» , Wa>» , 
etc. until the local optimum is found, (vii) Finally, mea- 
sure the difference between each fitness of the neighbors 
and the local optimum fitness. This provides us a rugged 
fitness landscape of the system. 



The above procedure starts from the particular initial 
state with the set of the random numbers Wi. Since we 
can change either the initial state to a different state 
in the 2 N states or the set to a different set chosen 
randomly, we can generate many samples. Each sample 
results in a different rugged fitness landscape of the 
system. Hence, we obtain an ensemble of them. Thus 
we can study the statistics of the structures of rugged 
fitness landscapes In this paper we take a thousand 
samples for the purpose. 
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TABLE 1. The relationship between the (real) output 
Wi and the (binary) inputs, {a tl , a i2 , • • ■ , a lk . , a.,}. 
Since there are (ki + 1) <7,'s, each of which has or 1, 
there are 2 ki+1 ways of inputs. These provide 2 ki+1 uj i , s, 
each of which is a real number randomly drawn from the 
interval [0, 1], according to a homogeneous distribution 

To apply the above method to the SF Kauffman 
model, we have to specify a model for the scale-free net- 
work with an arbitrary degree of nodes (k) — n, even in- 
teger. For this purpose, let us adopt a slightly modified 
version of the so-called Albert-Barabasi (AB) model[l2T] 
as a prototype model. In our model, we initially start 
with mo = n/2 nodes for seeds of the system, all of which 
are linked to each other such that the total link number 
is n(n — 2)1%. And every time when we add one node to 
the system, m — n/2 new links are randomly chosen in 
the previously existing network, according to the prefer- 
ential attachment probability of IL(fcj) = fc,/ X^^i* kj. 
Then, after t steps, we obtain the total numbers of nodes 
N(t) =n/2 + t and of links L(t) = n(n - 2)/8 + (n/2)t, 
respectively. We continue this process until the sys- 
tem size N is achieved. Hence, by this we can define 
(k) = 2L(t)/N(t) = n as t — ► oo. The generalization can 
be straightforward. Now we apply the above-mentioned 
dynamics to this modified AB model and call the result 
the SF Kauffmann model 0. 



IV. NUMERICAL RESULTS 

A. Rugged fitness landscapes 

Fig.l shows the rugged fitness landscapes of both the 
'quenched' random NK and the 'quenched' SF Kauff- 
man models. Here we have calculated the systems up to 
N = 256. This is just because of our computer power at 
this moment. We only show the data for N = 256 in this 
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FIG. 1: (Color online) Rugged fitness landscapes of the 
'quenched' random NK and of the 'quenched' SF Kauff- 
man models, where the total number of nodes is given by 
N = 256. (a)-(d) are shown for the random networks of 
K — 2,4, 8, 16, respectively, and (e)-(h) for the scale-free 
networks of (fc) = 2, 4, 8, 16, respectively, where (k) means 
the average degree of nodes. In each figure, the vertical axis 
shows the fitness while the horizontal axis shows the family of 
all one-mutant variants n omv of a local optimum. The dotted 
line denotes the fitness level of the local optimum and the 
curve the fitness differences between the local optimum and 
all the one-mutant variants. 



paper. As previously noted by kaiili'nian .'> the fitness 
landscapes of the random networks are very rugged. We 
see that the fitness landscapes of scale-free networks are 
very rugged as well, but quite different from those of the 
random networks. This can be understood as follows: In 
the random network each one of the nodes always meets 
with the same K links to the neighbors. The ruggedness 
can be dominantly bounded by the value of K. Hence, 
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FIG. 2: (Color online) The number of times h when a par- 
ticular value of the fitness difference D between a local opti- 
mum and n omv occurs in the rugged fitness landscapes. Here 
the 'quenched' models have been adopted, (c) and (d) [(a) 
and (b)] show the histograms and the tails for the scale-free 
[random] networks of (k) = (K =)2 (O), 4 (□), 8 (A), re- 
spectively, where N = 256. Here (k) stands for the average 
degree of nodes. The data for 1000 samples are superimposed 
in each figure, (b) and (d) (inset) are shown in the log-log 
(semi-log) plot of the tails. 



as K becomes large, the fluctuations in the rugged fit- 
ness landscapes become large. On the other hand, in 
the scale-free network there are various kinds of degrees 
hi of nodes. In other words, each node has its own fa 
links and there is the distribution of the degrees such as 
P(k) oc fc -7 . The AB model exhibits power-law with 
7 = 3. Therefore, the fitness landscapes can fluctuate, 
according to the degree distribution of the network. So, 
as the degree fa of a node is large, the difference in the 
fitness between the optimum and the mutants is expected 
to be large. Hence, the fitness landscapes obey the nature 
of the scale-free network. 



B. Histograms and tails of rugged fitness 
landscapes 

How can we detect the differences in the rugged fit- 
ness landscapes between the random and the scale-free 
networks? To do so, denote by h the histogram and de- 
note by D the fitness difference between the local op- 
timum and the one- mutant variant. Then we draw the 
histograms of the rugged fitness landscapes in the nor- 
mal plot [(a) and (c)], in the log-log plot [(b) and (d)], 
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and in the semi-log plot (insets), respectively (Fig. 2). 
Comparing (a) with (c) in the numerical results, we find 
that the histograms of the fitness landscapes for the ran- 
dom networks behave like Gaussian distributions, which 
can be fitted by h oc e~^ D ~ l ' D ^ l v whe re we se t the 
peak value as (D), the variance as v — U ((AD) 2 ) with 
AD = D — (D) . We have found numerically that (D) = 
{0.943,2.62,5.85} x 1CT 3 ; v = {4.31,5.75,7.58} x 1(T 3 
for K = 2,4,8, respectively. On the other hand, the 
histograms for the scale-free networks behave like non- 
Gaussian distributions with a broad tail, which can be 
fitted by 



(D-(-D)) z 
^ 



h oc 



if (D - {D)f « v 2 

e~ a « if D is comparable to (D) 

if D is much larger than (D) 



(4) 



Here we have found numerically that (D) — 
{1.48,2.74,5.27} x 1(T 3 ; v = {3.45,4.67,5.53} x 1(T 3 ; 
a = 1.30,1.21,0.978; = 3.11,3.50,2.98, for (k) = 
2, 4, 8, respectively. 

We note here the following: (a) The tail (i.e., scal- 
ing behavior) appears when the system size becomes as 
large as N = 256. And as the value of (k) is increasing, 
the value of (3 seems closer to 7 = 3 of P(k) in the AB 
model 12] [Fig. 2 (d)]. But in view of the limited accuracy 
of Fig. 2(d) (e.g., (3 = 3.50, 2.98), at the present moment 
this can be a conjecture, speculated from the numerical 
results for the system of N = 256. However, in Appendix 
A we give a sharp analytical arguments showing a strict 
relationship between the power-law decay and the "de- 
gree" fluctuations. And also, since we extend the system 
up to N — 1024, we are able to confirm ourselves that the 
tail behavior is maintained and become more prominent, 
as N is increasing. As an example, we show the result in 
Appendix B. 

(b) The results for the random networks show a kind 
of transition when the value of K is going up from K = 
2 to K = 8. This is consistent with the critical value 
of K c = 2 for p — 0.5 which was analytically obtained 
from the annealed model 5j. Therefore, the distribution 
below K c is quite different from that above K c so that 
the distributions for K > K c become more Gaussian- 
like as K increases. Very interestingly, we find a similar 
transition for the scale-free networks as well, when the 
value of (k) is going up from (k) = 2 to (k) = 8. 

This can be explained as follows: Suppose the dis- 
tribution of degrees is approximately given by P(k) = 
k~ 7 /C(l) such that we can impose normalization 
SfcLi •P(k) = 1, where £(7) is the Riemann's zeta func- 



tion defined by £(7) 
definition (k) — J2T= 



J2k=i h 1 ■ Substituting it to the 



J2k=i kP(k), then we obtain 
(k) = C( 7 - 1)/C(7), 



we obtain (k) = C(2)/C(3) = 1.64493 ■ ■■]!, 20205 •• • w 
1.3684. As studied by Aldana et. al.|8l Il9|. the critical 
value (k) c for the annealed dynamics of the SF Kauffman 
model is given by (k) c — C(lc~ 1)/C(7c) = 2 for p = 0.5 as 
well, where j c « 2.47875. Within the limited accuracy of 
our numerical results (some 20 percent) this is the same 
value as just stated for the quenched dynamics. In fact, 
we expect that also for our quenched dynamics there is a 
critical point around (k) c = 2 , with around 3, maybe 
again exactly at these values. In fact, in appendix A we 
show that also in our case the statistics of the rugged 
fitness landscape is bounded by K, such that (3 reflects 
the fluctuations of the 'degree', which should be the same 
both for the quenched and the annealed dynamics [8llT^|. 
Thus, our numerical results support this analytical result 
although our system is not very large but a finite scale- 
free network of N = 256. We also note here that we 
confirm that the phase transition at the critical value of 
(k) = {k) c = 2, which was proved analytically in the 
'annealed' model by Aldana et.al. 19], occurs in the 
'quenched' model as well p3| . 



C. The correlation between the fitness of local 
optima and their Hamming distance from the 
highest optimum found 
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FIG. 3: (Color online) The correlation between the fitness 
of local optima and their Hamming distance from the highest 
optimum found. It is shown for the random networks with 
K = 2 (a), 4 (b), 8 (c) and for the scale- free networks with 
(k) = 2 (d), 4 (e), 8 (f), respectively, where N = 256. The 
vertical axis stands for the fitness and the horizontal axis the 
Hamming distances between the largest local optimum and 
the local optima. The data for 1000 samples are superim- 
posed in each figure. In both cases there seem to exist phase 
transitions, i.e., at K — K c = 2 for the random networks and 
(k) = (k) c — 2 for the scale-free networks, supporting the ex- 
act statement for the 'annealed' version of the first mentioned 
case and the analytical prediction of Aldana et. al. on the 
annealed version of the SF model. 



which is finite for 7 > 2 and infinite for 1 < 7 < 2 
and which was first obtained by Aldana et. al.[H|. For 
example, since 7 = 3 for the special case of the AB model, 



Finally we present the correlation between the fitness 
of local optima and their Hamming distance from the 
highest optimum found 3] (Fig. 3). In both the random 
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and the scale-free networks we find the following: If K 
and (k) are as small as the critical value of 2, then the 
highest optima are nearest to one another. And the 
optima at successively greater Hamming distances from 
the highest optimum are successively less fit. Therefore, 
there is a global order to the landscape. On the other 
hand, as K and (k) increase, the correlations fall away. 
This shows that the previous assertions are maintained 
in the correlations, respectively. 



V. CONCLUSIONS 

In conclusion, we have studied the structure and statis- 
tics of the rugged fitness landscapes for the quenched SF 
Kauffman models, comparing with that for the quenched 
random NK Kauffman models. We have numerically cal- 
culated the rugged fitness landscapes, the distributions, 
the tails, and the fitness correlations of local optima with 
the Hamming distance from the highest optimum, respec- 
tively. From the results, we have concluded that in the 
SF Kauffman models there is a transition of network 
when (k) = (k) c = 2, while in the NK Kauffman mod- 
els such a transition occurs at K = K c = 2. This is, 
in some sense, quite analogous to the situation in the 
study of Boolean dynamics of NK and SF Kauffman 
models |8|- It would be very interesting if we could apply 
this approach to study fitness landscapes of other net- 
work systems. 
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APPENDIX A: ANALYTICAL ARGUMENTS OF 

RELATION BETWEEN THE POWER-LAW 
DECAY AND THE 'DEGREE' FLUCTUATIONS 

The tail behavior of the histogram is interpreted as 
follows. 

Define a network with N nodes, which can be any kind 
of network such as random, scale-free and exponential- 
fluctuating networks, where the i-th node has ki degree 
(i.e., link number). Let us consider the Kauffman model. 
As given in Table 1, the number of inputs to the i-th node 
is given by ki + 1. Let us now choose an iV-dimensional 
initial state A such as (1, 0, 0, 1, 1, • • • ,1). Considering 
the input from the given network structure, we can define 
fitness Wi on each node. Therefore, we can define the 
fitness Wa of the state A by 



W A 



1 N 

— yr 

N 2—* i 



(Al) 



where means fitness at the j-th node in the state 
A. Consider one-mutant family of the initial state A, in 
which there are N neighbor states with Hamming dis- 
tance one. For example, let us say one of them, B and 
define as (0, 0, 0, 1, 1, • • • , 1). In this example, only the 
state of node 1 is different from 1 to 0. Therefore, the 
difference in input values between the initial state A and 
this state B comes from nodes linked to node 1. This 
situation provides a difference in fitness. 

Suppose that the degree of node 1 is k\ and denote the 
ki nodes linked to node 1 by j% , ■ ■ ■ , ■ The values of 
nodes linked to node 1 are given different random num- 
bers W[ according to Table 1. We then obtain fitness for 
the state B as 



1 - i 



+ ■•■)• 

(A2) 

Thus, a genetic mutation in the state of one-mutant gives 
a change only for the node that the mutation occurred 
and the nodes linked to it. Therefore, if we consider only 
the fitness difference from the local optimum, then the 
fitness value of the one-mutant family that has Hamming 
distance 1 from the local optimum state depends upon 
which node the mutation occurs. Hence, in the case that 
there is mutation on node m, we obtain 

AW B (m) = i(AWi + AWji + • • • + AW 3hi ), (A3) 



where AW~ 



N 
W 

3 



Wj. The left hand side of Eq.(A3) 



means the fitness difference, D. To see what it means, 
let us define the averaged fitness difference for node m: 



(D B {m)) = 



We then have 



f 



-(AWt + AWj. 



AW jk ). (A4) 



D m =AW B (m) = -(D B (m)){k m + l). (A5) 

From this, if the average (D B (m)) is constant, then D m 
is proportional to k m + 1. But, more generally, there 
are two contributions: one from random number Wi in 
AWi+AWj 1 +■ ■ - + AWj km and another from degree k m . 
Here, if the random numbers are defined by a uniform dis- 
tribution, then we can understand that they contribute 
to the exponent of the fitness distribution function and 
the tail of the distribution function comes from that of 
the degree (link number). Because, since the maximum 
value of Wi is 1, it is bounded as 



|AWi + AW h 4 
Hence, this provides 



AW jhm | < k m + 1. (A6) 



\D m \ < (A7) 



N 

From the above, in the NK model the statistics of rugged 
fitness landscapes is bounded by K. In the SF model, 
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since ki is distributed by a power law, the statistics be- 
comes the same power distribution. Similarly, in the ex- 
ponential fluctuation distribution, so is the fitness dis- 
tribution. In this way, the statistics of fitness is strongly 
dominated by that of the link distribution in the network. 

APPENDIX B: TAIL BEHAVIOR OF THE 
SYSTEM OF N = 1024 

We show the tail behavior of the histogram in the sys- 
tem size of N = 1024 in Fig. 4. This may support our as- 
sertion in the text. We find again 7 rs 3 ; but since here 
the power-law decay is already obtained for < k >= 2, 
we cannot exclude that< k > c < 2, although such a state- 
ment would perhaps only reflect finite-size effects. 



FIG. 4: (Color online) log-log plot of the tail behavior 
of the histogram for scale-free networks of the system size 
N — 1024, where (k) = 2. Here the 'quenched' models have 
been adopted. The vertical axis means log of the number of 
times h when a particular value of the fitness difference D 
between a local optimum and n omv occurs in the rugged fit- 
ness landscapes. The horizontal axis means log of the fitness 
difference D. The data for 100 samples are superimposed in 
figure. We confirm ourselves that 7 » 3. 
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